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ABSTRACT 

We describe three-dimensional general relativistic magnetohydrodynamic simulations of a geometrically 
thin accretion disk around a non-spinning black hole. The disk has a thickness h/r ~ 0.05-0.1 over the radial 
range (2-20)GM/c 2 . In steady state, the specific angular momentum p rofile of the inflowing mag netized gas 
deviates by less than 2% from that of the standard thin disk model of Novik ov"& Thornel d 19731) . Also, the 
magnetic torque at the radius of the innermost stable circular orbit (ISCO) is only ~ 2% of the inward flux of 
angular momentum at this radius. Both results indicate that magnetic coupling across the ISCO is relatively 
unimportant for geometrically thin disks. 

Subject headings: X-ray: stars — binaries: close — accretion, accretion disks — black hole physics 



1. INTRODUCTION 

The recent development of general relativistic magneto- 
hydrodynamic (GRM HD) codes (e.g., iGammie et alJ 120031 
iDe Villiers & Hawlevll2003l) has finally allowed realistic nu- 
merical simulations of magnetized accretion disks around 
black holes (BHs). This has led to a bett er understanding 
of the inner regions of these disks (e.g. iKrolik et alJl2005[) 
and o f their role in launching relativistic jets (e.g. Mc Kinnevi 
2006). One of the interesting new results is the recognition 
that magnetic fields alter the structure of the accretion flow 
near and inside the innermost stable circular orbit (ISCO). 
This may result in large deviations fro m the tradit i onal pic- 
ture of a vanishing torqu e at the ISCO (Krolik 1999; Gammie 
[19991: IKrolik et alJIIOOl 

iPaczvnskil d2000h and lAfshordi & Paczvnskil d2003l) sug- 
gested that the zero-torque condition is likely to be a good ap- 
pr oximation for geom etrically thin disks. This was confirmed 
bv lShafee etafl d2008l) who, using a global height-integrated 
model, showed that modifications to the stress profile are neg- 
ligibly small for disk thicknesses h less than about a tenth 
of the local radius r. Their work was, however, based on a 
hydrodynamic model with a-viscosity and did not explicitly 
include magnetic fields. 

Numerical MHD simulations by [Krolik & Hawlevl d2002l) 
using a pseudo-Newtonian potential, and by IKrolik et alj 
(2005) using a GRMHD code, indicated that magnetic torques 
are indeed important inside the ISCO. In particular, these au- 
thors found no evidence for a "stress edge," leading them 
to argue that si mple disk models based on the zero-torque 
condi tion (e.g.,[Shakura & Sunvaevll973l : lNovikov & Thornel 
[19731 hereafter NT73) may be seriously wrong. If true this 
would undermine recent efforts to estimat e of the spin pa- 
rameters of BHs using the NT73 model dShafee et alj|2006t 
iMcClintock et al.ll2006b iDavis et alJl2006llLiu et alJl2008l) . 

The MHD simulations carried out so far have consid- 
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ered non-radiating accretion flows that are geometrically 
rather thick. The one exception is the recent work of 
[Reynolds & Fabian (2008) who considered a thin disk with 
h/r ~ 0.05 in a pseudo-Newtonian potential. In this paper, 
we describe global 3D GRMHD simulations of a thin disk 
(h /r ~ 0.05) around a non-spinning BH and compare our sim- 
ulated model with the NT73 model. 

2. NUMERICAL MODEL 

We use units with G = c= 1, e.g. the horizon is at r = 2M and 
the ISCO is at 6M. We report results in Boyer-Lindquist (BL) 
coordinates (f , r, 8, </)), referred to as the coordinate frame. In 
the expressions below g refers to the determinant of the met- 
ric. The fluid 4- v elocity and the magnetic field 4-vector (see, 
e.g.. lAnild[T989l for definitions) are given by m m and b^, and 
the rest-mass density, internal energy density, thermal pres- 
sure, and magnetic pressure as measured in the fluid comov- 
ing frame are p, u g , p = (7- \)u g with adiabatic index 7 = 4/3, 
and pb = b 2 /2. The total pressure is p tot = p + pt,. 

Simulati ons were performe d using a 3D GRMHD code 
HARM dGammie et all 120031) in Kerr-Schild coordinates 
using the interpolation scheme de scribed by iMcKinnevI 
(2006), inversion schem e described by lNoble et alJ (2006) and 
Mignone & McKinnev (2007), and other advances described 
by iTchekhovskoy et all d2007l) . For our fiducial model we 
used a 3D grid with resolution 512 x 128 x 32 corresponding 
to: (i) 512 cells in r, logarithmically spaced from r = 1.8M 
to 50M with "outflow" boundary conditions; (ii) 128 cells in 
8 going from to ir, non-uniformly spaced so that roughly 
half the cells are concentrated in the disk 5 ; (iii) 32 cells 
in <j), uniformly spaced from to ir/4. Cells at the disk 
equator had physical sizes roughly in the ratio of 2:1:7 in 
dr, rdd, rsm(8)d<j), which ensured that the turbulence was 
roughly isotropic in Cartesian coordinates, optimal for accu- 
rately resolving the magnetorotational instability (MRI). To 
test convergence we also used resolutions of 256 x 64 x 16, 
256 x 64 x 32, and 256 x 64 x 64. We also earned out 2D sim- 
ulations with resolution of up to 2048 x 256 x 1, but we do not 
discuss these results because turbulence decayed on an orbital 

5 We used a grid given by 0(jt (2) ) = [h(2x (2) - 1) + (1 -/i)(2.v <2) - 1) 7 + 1]tt/2 
for code coordinate x (2 ' with h = 0.15, giving roughly 6 X more angular reso - 
lution compared to equation (8) with h = 0.3 in McKinney & Gammie 1 2004). 



2 



time-scale (i.e. Cowling's anti-dynamo theorem holds). 

We began the simulation with an equilib rium torus 
dChakrabartil 119851; iDe Villiers & Hawlevl [2003b with inner 
edge at r = 2QM and pressure maximum at r = 35M and ad- 
justed the model parameters so the torus had h/r = 0.1 at 35M. 
We found that placing the torus at a smaller radius led to 
results too sensitive to the initial mass distribution. We de- 
fine h/r to be the density-weighted root mean square angular 
thickness of the disk at any given r, i.e. 
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(6). We also 



or written as A9„ 

consider the mean density- weighted thickness: h/r = A# a b s = 
(|A0|>. 

We embedded the torus with a weak magnetic field with 
j3 = p/pi, ~ 100. The initial field consisted of two poloidal 
loops centered at r = 28M and 38M to model a disorganized 
field with no net flux. (During the simulation, there is no 
organized flux threading the disk but some organized flux 
threads the BH.) The field strength was randomly perturbed 
by 50 % to seed the MRI instability. R ecent GRMHD simula- 
tions dMcKinnev & Naravanll2007allbl: iBeckwith et al.ll2008al) 
indicate that the results for the disk (but not the jet) should be 
roughly independent of the initial field geometry. The MRI is 
initially resolved with much of the torus having 10 cells per 
wavelength of the fastest growing MRI mode. 

In order to keep the accretion disk thin, an ad hoc cool- 
ing/heating function was added to the energy-momentum 
equations as a covariant source term (-(/' du g / 'dr) with 



dT 



7"cool 



where t is the fluid proper time. The gas cooling time (r coo i) 
was set to 2-k/SIk, where Qk = (r/M)~ 3 ' 2 M is the Keplerian 
frequency. Thus the gas was driven towards M eq , which we 
defined as that value of u g for which the specific entropy of 
the gas would be equal to the constant specific entropy (e.g., 
p/ p 1 ) of the initial solution. 

The simulation ran for a time of 10000M, corresponding to 
108 orbits at the ISCO (r ISC o = 6M) and 18 orbits at the initial 
inner edge of the torus. The results reported here correspond 
to averages computed over the period 6000M- 10000M, when 
the accretion flow had reached a quasi-steady state inside a 
radius of about 10M. Figure [TJ shows that the simulated disk 
had a thickness (given by eq. [T) of h/r ~ 0.06 -0.10 over the 
radius range of interest. The mean absolute thickness A6* a b s 
was smaller: ~ 0.04-0.07. At the ISCO, the two definitions 
of thickness gave h/r ~ 0.08, 0.06, respectively. 



3. RESULTS 

The flux of mass and specific angular momentum are given 



by 



M(r,t)-. 
L(r,t) 



pu 



(3) 



=4i(r,0-4ut0*,0, (5) 
respectively, where T£ is the r-<f> component of the stress- 
energy tensor, and 
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FIG. 1 . — Variation of the root mean square disk thickness h/r = A0 rms 
(eq. fj} as a function of radius r (solid line). Simulation results were averaged 
over the time interval 6000M— 10000M. Also shown is the mean absolute 
thickness A0 a (, s (dashed line). The vertical dotted line corresponds to the 
position of the ISCO. The initial torus had its inner edge at r m = 20M. 




T Lm = (p+u g + P + b 2 )u r u < p, 7^ out = b r b 4 



(6) 



FIG. 2. — Time-averaged profile of the specific angular momentum t m of 
the inflowing magnetized gas (solid line), compared to the idealized profile 
assumed in the NT73 disk model (dashed line)._The horizontal dotted curve 
shows the net flux of angular momentum t- m - £ ut ■ This is independent of r 
up to ~ 10M, indicating the simulation has achieved steady state well-beyond 
the ISCO. 

The specific energy flux E/M is found by replacing T£ with 

T t r in the above equations. In equation (O, the term l m rep- 
resents the specific angular momentum advected inward with 
the accretion flow. Similarly, £ out represents outflow of spe- 
cific angular momentum as a result of the magnetic shear 
stress. In steady state, each of these terms is independent of t, 
and their sum is independent of both r and t . 
Figure|2]compares the time-averaged profile of l- m from the 
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FIG. 3. — Lower panel shows time-averaged profile of the fluid shear stress 
(solid line) compared to the NT73 stress (dashed line). Upper panel shows a 
viscosity coefficient (solid line). The non-zero stress in the plunging region 
causes only a minor deviation from the NT73 result for the flux of specific 
angular momentum. 

simulation, averaged over an angular range SO = ±0.2 around 
the disk mid-plane (2 to 3 density scale heights), with that 
predicted by NT73. In the latter, £- m is equal to the Keplerian 
specific angular momentum for all radii down to the ISCO; 
inside the ISCO, £- m is taken to be constant since, by assump- 
tion, no angular momentum is removed from the gas in the 
plunging region. The £- m profile from the simulation shows 
modest deviations from this idealized_profile: (i) £ m is slightly 
sub-Keplerian outside the ISCO; (ii) continues to decrease 
for a range of r inside the ISCO; (iii) £- m becomes essentially 
independent of r close to the horizon. 

The surprising feature of Fig. 2 is that the two curves are 
so close to each other inside the ISCO. While it is true that 
the profile of £- m from the simulation drops inside the ISCO 
to a value of 3.39M, this value is only 2.0% less than the 
NT73 value of 3.464M. Thus, the magnetic coupling that 
technically operates inside risco does not have much real ef- 
fect on the accreting gas. At the ISCO £ oat = 0.075M, which is 
~ 2.2% of 4,. Further, E/M at the horizon is 0.940, which is 
only slightly different from the NT73 value of 0.9428, show- 
ing that the accretion efficiency is not enhanced by magnetic 
fields. For two-dimensional simulations of thicker disks with 
hjr ~ 0.3, iMcKinney & Gammid (12004 found that L/M = 
3.068 and E/M = 0.950. It appears that L/M deviates from 
the NT73 value roughly proportional to h/r and so thin- 
ner disks should be even closer to NT73. These results are 
converged since, between a resolution of 256 x 64 x 32 and 
512 x 128 x 32, the value of L/M changes by less than 0.2% 
and E/M changes by less than 0.5%. 

Figure [3] shows the time-average of the normalized shear 
stress 

W =10 3 J j ' Tfry/=gd0d<l>/{2nrM), (7) 

where W = 10" 3 WM as equation (5.6.1b) in NT73 and T^ f is 
the orthonormal stress-energy tensor components in the co- 
moving frame. We restrict the integral to fluid with u,(p+u + 



FIG. 4. — Time-averaged profile of £turb, which measures the relative mag- 
nitude of turbulent fluctuations in the accreting gas. The fluid becomes mostly 
laminar inside the ISCO. 



p + b 2 )/p > -1, i.e., only gravitationally bound fluid (no disk 
wind); this lowers the value of W near the horizon by 50%. 
Figure|3]also shows the space-time averaged viscosity param- 
eter a = Tiff that is roughly 0.1-0.15 outside the ISCO. 
Between the ISCO and r = 3M, a rises to 0.6 and finally drops 
to ~ near the horizon. There are non-vanishing stresses 
inside the ISCO and the "stress edge" is near the horizon. 
These stresses appear due to organized magnetic flux, but ev- 
idently do not cause significant transport of angular momen- 
tum (Fig. |2j and are not expected to be a significant source of 
dissipation and radiation. 

Another way of evaluating the degree of coupling is to con- 
sider, by analogy with steady state flows, the locations of 
"critical points," defined as the radii at which various volume- 
averaged outgoing characteristic speeds vanish in the coor- 
dinate frame. For our GRMHD simulation, the fast magne- 
tosonic radius is located at rf ast ~ 5.1M showing that, even for 
this relatively thin accretion disk, magnetic fields do enhance 
the communication between the disk and the plunging region. 
On the other hand, the Alfven critical point is much further out 
at 8.0M. Since Alfven waves play an important role in angu- 
lar momentum transport, this suggests that the shear c oupling 
via ma gne tic fields is wea k, unlike the suggestion of iKrolikl 
(119991) and lGammid (11999b . 

For fitting disk models to observations the most useful 
quantity is the rate of dissipation of energy as a function of 
radius, since this is what determines the radiation emitted by 
the accretion disk. As a simple estimate, consider the stan- 
dard theory of viscous disks in which the energy dissipation 
rate is equal to the local shear stress multiplied by rd^l^/dr 
(this is not necessarily valid for an MHD fluid, and is even 
less valid in the plunging region). In contrast to the standard 
model, which assumes a vanishing torque at the ISCO, we 
find £ oul = 0.022^i„ at r = nsco- This modification to the in- 
ner boundary condition will cause the luminosity of the disk 
outside the ISCO to increase by about 4%. 

Figure |4] shows one measure of the magnitude of turbulent 
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fluctuations, 



„r\2 



(8) 



time and angle-averaged as for Fig. [2] We see that the gas 
at radii > 10M is highly turbulent due to the MRI. However, 
the gas flow is fairly smooth at smaller radii, suggesting that 
there is little dissipation of turbulent kinetic energy inside the 
ISCO. 

4. CONCLUSION 

For a geometrically thin accretion disk with h/r ~ 0.05-0. 1 
(Fig. 1) around a non-spinning BH, we find that the specific 
angular momentum profile £ m {r) of the accreting magnetized 
gas is quite similar to that assumed in the idealized model of 
NT73 (Fig. 2). Specifically, l- m is slightly smaller than, but 
nearly equal to, the Keplerian value of £ for radii outside the 
ISCO, and is almost independent of r inside the ISCO. At the 
BH horizon, the value of l- m deviates from the NT73 value by 
only 2.0%. 

This result suggests that the NT73 model is a good approx- 
imation for thin disks. However, fast magnetosonic waves 
from radii as small as 5.1M can escape the plunging region 
and communicate magnetic stresses between the disk and the 
plunging regions, which means that the zero-torque condition 
assumed in the NT model is not perfectly valid. Nevertheless, 
the normalized outward flux of angular momentum £ out at the 
ISCO is only ~ 2.2% of the inward flux l- ln . Thus, the mag- 
netic coupling across the ISCO is quantitatively weak. As a 
result, we estimate the luminosity of the disk to be enhanced 
by no more than a few percent. 

We find that turbulent activity is pronounced at radii be- 
yond about 10M, and that the flow is nearly laminar inside 
the ISCO. This suggests that the bulk of the dissipation oc- 
curs outside the ISCO in the disk proper, not in the plung- 
ing region. In our next paper we plan to discuss the energy 



dissipation profile of 3D GRMHD simulations, which is ulti- 
mately what determines the observed radiation. Future studies 
should also investigate the dependence of these results on hj r, 
magn etic field geometry, black hole spin (e.g. Gam rnie"et alj 
2004), cooling, resistivity and viscosity, and mass distribu- 
tion. 

The results reporte d here are in qualit a tive agreement with 
those obtained by iRevnolds & Fabian] d2008l) They car- 
ried out a 3D non-relativistic MHD simulation in a pseudo- 
Newtonian potential, whereas ours is a fully relativistic 
GRMHD simulation. Also, their initial torus had its inner 
radius at the ISCO, r m = 6M, whereas our torus initially had 
rj n = 20M. This might explain some differences in our results, 
e.g., their turbulence edge is at 5.5M whereas ours is at 6.8M 
(compare their Fig. 3 with our Fig. 4). Indeed, we find that our 
results are sensitive to the initial conditions if we place the in- 
ner edge of the torus much inside r m ~ 20M. Despite these 
differences the two simulations agree on the main result, viz., 
for a geometrically thin disk, the ISCO does behave like a 
physical b oundary separating the d isk proper from the plung- 
ing region. Beckwith et al. (2008b) suggest a larger deviation 
from NT73, but they considered a torus with a larger value 
of h/r and did not use a self-consistent dissipation model. A 
study of the actual dissipation in the simulation, the goal of 
our next paper, will hopefully resolve these differences. 
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